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ABSTRACT 

We present a ray-tracing technique for radiative transfer modeUng of complex 
three-dimensional (3D) structures which include dense regions of high optical 
depth like in dense molecular clouds, circumstellar disks, envelopes of evolved 
stars, and dust tori around active galactic nuclei. The corresponding continuum 
radiative transfer problem is described and the numerical requirements for in- 
verse 3D density and temperature modeling are defined. We introduce a relative 
intensity and transform the radiative transfer equation along the rays to solve 
machine precision problems and to relax strong gradients in the source term. For 
the optically thick regions where common ray-tracers are forced to perform small 
trace steps, we give two criteria for making use of a simple approximative solver 
crossing the optically thick region quickly. Using an example of a density struc- 
ture with optical depth changes of 6 orders of magnitude and sharp temperature 



- 2 - 



variations, we demonstrate the accuracy of the proposed scheme using a com- 
mon 5th-order Runge-Kutta ray-tracer with adaptive step size control. In our 
test case, the gain in computational speed is about a factor of 870. The method 
is applied to calculate the temperature distribution within a massive molecular 
cloud core for different boundary conditions for the radiation field. 

Subject headings: radiative transfer — methods: numerical — stars: formation 
— (ISM:) dust, extinction — ISM: globules — galaxies: active 

1. Introduction 

Almost all our information about astrophysical objects have been obtained by analyzing 
their radiation. Therefore, radiative transfer (RT) modeling of spectra and images is a 
standard technique of astrophysical research. Moreover, radiation transport is an important 
physical process to transport energy and momentum, often controlling the energy balance 
of the object, or altering its appearance by radiation pressure (see, e.g., Yorke & Sonnhalter 
2002; Richling & Yorke 2000). Given those two major applications of radiative transfer, 
it could be expected that radiative transfer modeling is a well-advanced and elaborated 
astrophysical tool being established for many years already. On the contrary, the solution 
of 3D radiative transfer problems is still one of the outstanding challenges in computational 
astrophysics. This is due to the high dimensionality of the radiation field (coordinates for 
the location in space, for the direction, the wavelength, and the time), the many orders of 
magnitude in the change of density and spatial scale that may have to be covered, the integro- 
differential character of the transport equation, and the nonlinear coupling of the different 
energies ranges through the local energy balance (see, e.g. Steinacker et al. 2003). Together 
with the many free parameters entering complex structures, this also has prohibited modeling 
of images with inverse radiative transfer, aside from very simple cases where semi- analytical 
solutions can be found (Steinacker et al. 2002). 

The currently proposed methods to solve the continuum RT problem are based on the 
Monte-Carlo approach, finite differencing schemes on (adaptive) grids, moment approaches, 
and solving the radiative transfer equation along rays (ray-tracing) , often with a mixture of 
techniques implemented into one program (for a review of radiative transfer techniques, see 
e.g. Pascucci et al. 2004; Henning 2001; Kudritzki & Hummer 1990; Hubeny 1990; Mihalas 
1978). 

A common problem of most RT codes is the treatment of the densest regions of as- 
trophysical objects like molecular cloud cores, accretion disks, or dust tori around active 
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galactic nuclei, where the optical depths can reach high values like 10^ depending on the 
wavelengths of the considered radiation (see Tab. 1 for typical peak values). Applied to 
these objects, Monte Carlo simulations have to calculate many interactions of the photons 
(for an improvement see Wolf 2003), adaptive grid techniques refine the regions into many 
cells (Steinacker et al. 2002), and ray-tracer are forced to do small tracing steps. Moment 
methods are well-posed to treat the optically thick regime. However, in parts of the com- 
putational domain where the optical depth is at the order of or smaller than unity, the 
radiation field can become strongly varying with direction, and many orders of momenta are 
required to describe the radiation field correctly. Moreover, errors in calculating the optically 
thick regions correctly have little infiuence on the energy balance within the computational 
domain, but will alter the appearance of the object strongly as the outgoing radiation of 
the optically thick regions defines the inner boundary for the layer of the object where the 
optical depth is of the order of unity. 

Ray-tracing solvers have the advantage that the numerical error can be controlled pre- 
cisely. General purpose schemes to solve ordinary differential equations along rays are avail- 
able in high-order accuracy and with adaptive step size control, like the advanced 5th-order 
Runge-Kutta solver proposed in Press et al. (1978). The drawback of such solvers, however, 
is that the known functions of the transport equations have to be evaluated very often in 
order to cancel all low-order error terms. Aside from this, they are not making exphcit use of 
the fact that the radiation field becomes very simple when the optically thick case is reached. 

The ray-tracing follows the transport of the radiation often given on an underlying grid. 
In the short characteristics approach, the change of the intensity is interpolated onto the 
local grid (see, e.g., Mellema et al. 2006), while in long characteristics approaches, radiation 
is transported between two points. 

For the case of the spectral energy distribution of a circumstcllar disk around a low- 
mass star, different numerical RT techniques have been compared in Pascucci et al. (2004) to 
estabhsh a benchmark for 2D continuum RT. A somewhat unreahstically fiat radial profile 
with a power law index of -1 for the density distribution and a large inner radius was assumed 
to reduce the optical depth in the disk and therefore the numerical effort. Steepening the 
radial profile towards more realistic values around -2 and assuming inner disk radii around 
a few stellar radii up to 10 stellar radii will increase the optical depth in the midplane to 
values around 10^. The next step, therefore, is to develop numerical schemes that are able 
to perform accurate calculation even in region where the optical depth is larger than 100 
(maximal value in the 2D RT benchmark Pascucci et al. 2004). In this paper, we suggest a 
scheme for ray-tracing to overcome this problem. 

We propose to transform the ray equation into a form that is more appropriate for 
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high-order ray-tracers applied to high-opacity structures. In Sect. 2, we define the radiation 
transport equation and the intensity transformation as well as the modified transport equa- 
tion. We discuss the advantages of the transformation in terms of limited machine precision 
and applicability of approximative tracing in the optically thick regions. We apply the solver 
to the test case of radiation crossing a medium with large density and temperature gradients 
in Sect. 3 and compare the results to the solution obtained with common ray-tracing. In 
Sect. 4, we use the density structure from a recently performed multi- wavelength modeling of 
a complex dense molecular cloud core, rescaled to higher masses, to investigate the thermal 
infiuence of a nearby illuminating star. The findings are summarized in Sect. 5. 

2. Transformation of the continuum radiative transfer equation 

While the scheme proposed in this paper is applicable to both line and continuum RT, 
we will introduce and demonstrate its capabilities by concentrating on continuum radiation. 
We consider a configuration of gas and dust being exposed to radiation from internal or 
external sources (for astrophysical applications, these might be stars, nearby gas and dust 
etc.). The radiation may interact with the dust particles in form of absorption or scattering, 
and the dust particles are treated like black body emitters with the temperature T. The 
stationary specific intensity Ix of the radiation field is a function of the location x, the 
direction n, and the wavelength A. The corresponding RT equation reads 



with the phase function p describing the probability to scatter radiation from direction n' 
into n and the specific Planck function Bx at temperature T. The absorption and scattering 
coefficients k"^* and k*"^" are defined by the product of the number density at the location x 
and the cross section for absorption and scattering, respectively. Scattering of radiation at 
dust particles is important for short wavelengths and depends on the size, shape and chemical 
composition of the dust particles. The scattering integral in (1) makes the application of an 
explicit ray-tracing scheme impossible, and iterative schemes are commonly used by fixing 
the intensity to a known start value Iq, calculating the scattering integral, and solving (1) 
for an updated Ji until convergence is reached. Another way is to split the intensity into 
an unprocessed, a scattered, and re-emitted component and to treat these components with 
different solvers (Steinacker et al. 2003). We will assume in this paper that the scattered part 
of the radiation is treated separately and will concentrate on the unprocessed and re-emitted 
radiation causing difficulties when regions of high optical depth are encountered. 



nVgIx{X,x,n) 



[n'^'^iX, x) + k'^{X, x)] Ix{X, f , n) + «"^^(A, x)Bx[X, T{x)] 




(1) 
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We trace the radiation along a ray with a given initial value Ix{X,xo) at the point xq 
where a; is the coordinate along the ray. The solution of (1) in the direction n of the ray for 
vanishing scattering is 

X 

h{X, x) = h{X, 0) exp [-t(A, 0, x)] + J dx' k'^^'^X, x')Bx[X, T{x')] exp [-t(A, x' , x)] (2) 



with the optical depth 

Xj 

T{X,Xi,xj) = [ dx' k'''''{X,x'). (3) 



Introducing a typical length scale L for changes in the intensity or density structure to 
be considered, wc define a region within the computational domain to have a high optical 
depth if k"'^^{X,x)L ^ 1. For the case of local thermal equilibrium, Eq. (2) approaches the 
limiting case I\{X,x) = Bx[X,T{x)] as the local thermal emission dominates the intensity 
due to the strong damping of emission at preceding locations along the ray. 

Although the differential equation is of simple structure and the solution incorporates 
only an integration over a known function, the numerical evaluation of (2) is not straight- 
forward. The reasons are the two exponential functions entering the radiation equation. 
First, the exponential containing the optical depth has to be calculated precisely along the 
ray. Given a limited computational range of the computer, an optical depth of 1000 usually 
exceeds these limits and makes it necessary to renormalize the intensity. Second, the Planck 
function 

rises sharply with A for temperatures in the Wien part T{x) < t(A), with t{X) = he/ Xk and 
6(A) = 2hc^/X^. 

A substitution of the intensity transforming the quantities into a range where numerical 
limits of the computer are not encountered and the exponential terms of the equation can 
be controlled is therefore desired. 

Another requirement for the ray-tracing is speed. To model complex density structures, 
a parameter number of several hundreds is usually required (sec, e.g., Stcinackcr ct al. 2005). 
Given a certain parameter set describing the density structure, the temperature within the 
object is determined and used to calculate theoretical images at different wavelengths. These 
images are compared with the observed images and a minimization tool like simulated an- 
nealing is used to optimize the model parameters. Assuming several thousand pixels to be 
fitted simultaneously, and a minimum of 10^ steps of the optimization tool, the total number 
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of rays often exceeds 10^. To perform the calculation within a reasonable time, the ray-tracer 
should be optimized with respect to speed. 

An advanced ray-tracer with 5th-order accuracy has been proposed by Press et al. 
(1978), based on the Runge-Kutta scheme. It is coupled with an adaptive step size control 
using an embedded form of the 4th-order Runge-Kutta formula. As parameters for the error 
truncation, values determined by Cash & Karp (1990) are used. 

The disadvantage of this tracer is that the density and the temperature function are 
called six times each in every tracer step. Each of the called functions contains several 
hundred parameters for a complex 3D structure making an evaluation expensive. On the 
other side, the tracer steps are chosen adaptively and with full error control. The tracer is 
therefore the first choice for astrophysical problems with moderate optical depths variations. 

A major numerical difficulty arises for regions with high optical depth «;"^*(A, x)L > 100. 
Neglecting scattering, the change in the intensity expressed in the right-hand side of (1) 
is proportional to /t"**(A,a:). Having calculated the intensity at a point x, the difference 
\Bx{X,x) — Ix{X,x)\ is small as the local thermal emission dominates the radiation field in 
the optically thick regime. A small change in the intensity suggests that the solver is able to 
perform large spatial steps. However, the difference \B\{X,x) — /a(A, a;)| is multiplied with 
the large absorption term k"'"*(A, x) amplifying any change in the source term that is arising 
from the spatial variation in the temperature. To meet the accuracy requirements, the tracer 
therefore performs small steps although the change in the intensity is known in the optically 
thick limit and the local source term has no strong gradients. A modification of the tracer, 
therefore, woiild be desirable in order to both crossing the regions of high optical depths 
quickly without introducing large errors and to keep the solution away from the fioating 
point limits. 

In order to develop a ray-tracer meeting those requirements and being acceptable in 
terms of computational effort, we introduce the relative intensity 

For a vanishing source function, it approaches unity, and for the optically thick part with 
Ix{X,x) !v Bx{X,x), Dx{X,x) has a value close to 1/2. Inserted into (1) and neglecting 
scattering, we obtain the differential equation 

f HI - ^.(A..)) (.-(A..) [. - 2.,,.)] - Ifl^) . (6) 

The relative intensity is chosen in a way that the source function in the differential equa- 
tion with the exponential dominating the Wien approximation regime is transformed to the 
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function 

1 dBx{X,x) _ t{X) dT 1 
BxiX,x) di ~ T{xyd^l-e-'W/^' 

The denominator 1 — exp[— t(A)/T] has a value close to 1 for t{X) 3> T so that the exponential 
terms no longer cause numerical problems. 

Another question is if the used rational function transformation amplifies numerical 
errors. However, no strong solution error amplification is introduced by transforming to the 
relative intensity, as the relative error in D\[X,x) 

±ADx{X,x) ^ ±AIx{X,x) 2hiX,x) + Bx{X,x) ^ ±2AhiX,x) 
Dx{X, x) h{X, x) h{X, x) T A/a(A, x) + Bx{X, x) " h{X, x) t A/a(A, x) ^ ^ 

is at most only about a factor 2 larger than the error in the intensity if A/;^(A, x) <^ h{K x). 

To obtain a criterion for the use of an approximate solver in the optically thick region 
Dx{X,x) ~ 1/2, we introduce a positive limit e for the relative difference of intensity and 
Planck function with 

\Bx(X,x)-Ix(X,x)\ \l-2Dx{X,x)\ 



Bx(X,x) l-Dx(X,x) 

Using Ist-order finite differencing for (1) for two neighboring points Xi and X2 and neglecting 
scattering, the limiting case of condition (9) reads 

^ ""ff, , = ±£ (10) 

AxK,'''''{X,x) ^ ' 

with Ax = X2 — X\ and ABx{X,x) = Bx{X,X2) — Bx{X,xi). The sign in front of e refers to 
the sign of Bx{X,x) — Ix{X,x). The condition is fulfilled either due to small relative changes 
in the source term |A_Ba(A, x)\/Bx{X, x) or a large local optical depth Axk,"''^'^{X, x). In a pre- 
calculation along the ray, the regions where the condition is fulfilled can be identified. Then, 
a fast solution can be obtained for these region by performing large steps while assuming 
Da(A, x) — 1/2. As the hmit for validity, we have used e — 10~^ in this paper. 



3. Application of the ray-tracing scheme to a configuration with 1 < r < 10^ 

To test the proposed ray-tracer, we choose a configuration with variations in the local 
optical depth covering six orders of magnitude. The temperature profile is chosen so that it 
includes both external heating and a deeply embedded radiation source. Such a configuration 
is often found, e.g., in star formation regions where young stars illuminate clouds where others 
embedded star are currently forming. 
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As density profile, a sum over two Gaussian distributions is used for the number density 
of the absorbing material (representing, e.g., two molecular clouds containing gas and dust) 



n{z) — rii exp 



z — Zi 



+ n2 exp 



Z - Z2 
W2 



(11) 



with the two optical depths (i = 1,2) Tj(A, x) = Ki{\,x)L — a'^^{\)ni{z)L having the 
values 10^ and 10"^, respectively, the absorption cross section (t"^*(A), and introducing the 
dimensionless spatial variable z = x/ L. The widths of the Gaussians are chosen as Wi = 1/2 
and W2 = 1- We will calculate the intensity along a representative ray that crosses both 
density maxima and encounters strong variations in the source term. In a real application, 
of course, all rays will be calculated (see Sect. 4). 

Without calculating the correct temperature reached in local thermal equilibrium for 
this test case, we assume a temperature distribution close to what can be expected for the 
density distribution defined in (11) 



T{z) = To{z + zo)-^'^ + T:,ex.v 



^3 



W2, 



+ T4 exp 



W4 



bg- 



(12) 



The first term represents a temperature profile arising from an external heating by a distant 
radiation source (e.g. a star) located in an optically thin region at large negative z (Tq = 50 
K, zq = L/2). Numerically, this corresponds to an exponential drop of the intensity that 
would quickly reach the machine precision limit without a source term. Therefore, the source 
term contribution, although being low, is determining the intensity in the center of the clump 
and saves the solver from failing due to precision problems. 

An intermediate region is reached between the two clumps where external radiation can 
illuminate the material from other directions. Therefore, we have chosen the temperature to 
have a broad Gaussian profile with T3 = 10 K, z^ — 3.7L, and — L, peaking where the 
density has a minimum. 

The second clump is less dense, but hosts an embedded source of heat represented 
by the second compact Gaussian profile in (12) with T4 — 1500 K, Z4 — 6L, and — 
L/10. In this case, the tracer enters a region where the source distribution rises sharply 
and dominates the radiation field again. A constant background radiation field in the far- 
infrared wavelength region will heat up the matter slightly, represented in our test case by 
the constant background temperature Thg = 5 K. 

We have chosen a wavelength of A = 6 //m for this test case because the thermal 
emission defined by the temperature profile defined in (12) gives a substantial contribution 
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to the radiation field. Note that in common star formation regions, the local optical depth 
will reach 10^ only in the densest regions at this wavelength. Fig. 1 gives an overview over 
the test case and the numerical solution of the corresponding radiative transfer equation. 
The upper left panel shows the local optical depth La"-^^ {X)n{z) as a function of z with the 
two Gaussian profiles establishing a density change of up to six orders of magnitude. In the 
lower left panel, the temperature variation along the ray is given, with a slow decrease in the 
outer parts of the first density clump, a local maximum between the clumps, and a strong 
bump in the center of the second density maximum. 

Applying the Runge-Kutta scheme with an solution accuracy of 10~^ and ~213000 main 
tracing steps to the entire region, we show in the upper right panel the numerical solution of 
the radiative transport equation /{""(A, z)/I\{\ 0) along the ray. The normalized intensity 
decreases exponentially from the start value 1, entering a region where the thermal emission 
dominates the solution at 2; ~ 0.4. The shape of the intensity reproduces the temperature side 
maximum at 2; = 3.8. At z ~ 4.8, the Planck function rises quickly due to the thermal peak 
centered at z = 6. The Gaussian shape in the intensity is broader than the corresponding 
temperature profile as the optical depth is r < 10'^ so that a small contribution of preceding 
radiation is added to the local thermal emission in the source integral of (2). An additional 
factor for the broadening is that for low temperatures T[z) < t, the Planck function rises 
as ~ exp[—t/T{z)], while around the peak ai z — 6 with t KiT{z) it is turning to follow the 
Bx(X, z) ~ T(z)-dependence for the large temperatures. 

The solution of the transport equation using the proposed tracer for the relative intensity 
Dx{\ z), transformed back to an intensity /j^'^'(A, z) using (5), is overlayed in the upper right 
panel and can not be distinguished from ll^'\X, z). 

The absolute value of the difference between the full solution and the relative solution 

Dxi\,z)Bxi\,z) 



(13) 



l-Dx{X,z) 

is shown in the lower right panel, normalized to the start intensity (thin line). The dotted 
vertical lines a.t z ^ 0.8 and z ~ 3.2 indicate where the relative tracer has switched to an 
accelerated crossing of the optically region and back to full tracing, respectively. 

Before discussing the solution error, it should be considered that calculating the differ- 
ence of the solutions requires to interpolate one solution onto the other, as they have been 
obtained using different adaptive step sizes. This introduces a new interpolation error into 
the difference that has to be distinguished from the error related to the solution process. To 
estimate the order of magnitude interpolation error, we have overplotted in the lower right 
panel the absolute value of the error of a linear interpolation of the full normalized intensity 
onto the z- values of the relative intensity values as a thick line. The second derivative of 
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Ix{X,z) can not be used, as the full solution is only known at discrete values, so that calcu- 
lating the second derivative would introduce new interpolation errors into the error analysis. 
We therefore estimate the error in the optically thick regions with Dx{X, 2;) ft! 1/2 to be 



2/a(A,0) dz^ 
with the second derivative of the Planck function being 



(14) 



dz'^ dz 



1 + e-*(^)/^(^) dBx{\ z) 2 dT{z) ^ £T{z) [ dT{z] 



Bx{X,z) dz T{z) dz dz^ \ dz 

(15) 

The interpolation error shown as a thick line in the lower right panel starts at a high value 
of 10~^ as the step size of the full tracer is large in the region where the density is low and 
the thermal contributions are small. The error becomes very small in the region of high 
optical depth as the full tracer is enforced to use small steps entering the interpolation error 
quadratically. The interpolation error has roots wherever the temperature has extrema as 
d'^Bx{X,z)/dz^ ^ dT{z)/dz. 

For small z, the plotted solution difference is of the order of the interpolation error, and 
the real difference is likely to be much smaller as both solutions are obtained by the 5th- 
order Runge-Kutta scheme in this region. In the optically thick region between the first two 
dotted lines, the interpolation error is below 10~^, and the thin curve shows the difference 
between the full scheme and the rapid cross stepping used by the relative solver. Still, this 
error is small and reveals the real benefit of the accelerated tracer. While the full solver 
spends 209000 grid points to trace the thick region, the relative solver uses 240 steps only 
with a difference in the solution that is of the order of 10""^ to 10~^. 

A second pair of dotted lines appears at ^ ~ 6, but the interval is so small that it 
appears as one line. The reason that the relative solver is not switching to larger steps 
although the local optical depth is 10^ for the second density maximum is related to the 
large gradient in the source function. Only in the maximum of the source term, the change 
in the source function becomes small enough so that the change in the relative intensity 
allows large stepping for a small range. The difference in the two full solutions is up to 
an order of magnitude larger than the 10~^-accuracy limit of the 5th-order Runge-Kutta 
scheme, as the error in D is twice the error in I, and interpolation errors and tracer errors 
are added up in the thin line. 

To summarize the findings for this particular test case, the relative tracer can reproduce 
the solution achieved by the full tracer with an accuracy of about 10~^ while it crosses 
the optically very thick regions 870 times faster than the full tracer. The relative tracer is 
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sensitive to changes in the source function to achieve this accuracy even in cases where strong 
intensity changes are caused by a varying source function within optically thick regions. 
The concrete gain in performance when using the relative tracer will depend on the extent 
and optical depth value of the dense regions to be crossed as well as on the occurrence of 
complications like embedded sources. 

We have performed 1000 test runs with varying total optical depth, density gradient, 
and source function gradient to explore the overall influence of the dominant ingredients 
on the run time and the error maximum. The resulting dependencies are shown in Fig. 2. 
The top right panel shows the relation between the total optical depth for a ray crossing the 
considered region and the run time normalized to the test case discussed in this section. They 
are directly proportional, as the adaptive steps of the Runge-Kutta solver scale inversely with 
the optical depth. For higher optical depth, the approximation of the relative tracer is better 
fulfilled, decreasing the maximal error for higher total optical depth. Therefore, the maximal 
error decreases inversely with the run time, as shown in the top left panel. The dependence 
on the density gradient is depicted in the lower left panel. For small run times, the run time 
rises linearly with maximal density gradient. This rise gets weaker for larger run times, as 
most of the run time is taken by the overall high optical depth. The lower right panel shows 
that there is no dependency of the run time on the gradient in the source function as the 
large optical depth dominates the run time in the region where the relative intensity solver 
is used. 



4. Application of the ray-tracing scheme to 3D illumination effects of a dense 

molecular cloud core 

To demonstrate the capabilities of the new ray-tracing scheme, we apply the scheme to a 
more realistic astrophysical problem, incorporating a large number of ray-tracing calculations 
through a high-opacity massive molecular cloud core. These cores are of particular interest 
as they are thought to be the progenitors of young stars. Characterizing their physical 
conditions will lead to a better understanding of how stars are forming in molecular clouds. 
The detailed understanding of images of cores crucially relies on the knowledge of the 3D 
structure, as projection effects limit the reliability of simple ID modeling (e.g., Steinacker 
et al. 2004). In Steinacker et al. (2005), we have fitted images of the cloud core p Oph D 
in the star formation region Rho Ophiuchi that have been obtained at mid-infrared (MIR) 
and mm-wavelengths. The complex elongated structure was modeled using a series of 30 3D 
Gaussian distributions, determining both the density and temperature structure by inverse 
multi-wavelength 3D continuum radiative transfer. 
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The density structure deviates significantly from spherical symmetry and has several 
local maxima. Hence, it is well-suited to serve both as a realistic test application for the 
ray-tracing scheme and as a representative shape of a dense core. The total core mass of 
2.3 solar masses, however, is not sufficient to reach substantial optical depth (>10) at MIR 
wavelengths. We therefore increase the core density uniformly by a factor of 100 until the 
optical depth ranges from 0.1 to 1000, placing the core among the most massive dense cores. 
The temperature structure for the high mass core is very different from the distribution 
determined for moderate optical depths and an outer radiation field of a low-mass star 
formation site like Rho Ophiuchi. It depends on the radiation field produced by nearby 
stars, which is most likely dominated by the massive stars in the considered star-forming 
region. The strong pressure of their emitted radiation and the stellar wind alter the structure 
of the parental molecular cloud substantially by opening up gaps and compressing the cloud 
material. 

In this section, we assume that the cloud core is still surrounded by a diffuse molecular 
cloud with an optical depth of 1 at A = 550 nm, and that a nearby B6 star at a distance 
of 1 pc is illuminating the core. The strong UV radiation of the star can not penetrate the 
parent cloud and is absorbed by gas, dust particles, and polyaromatic hydrocarbons (PAHs) 
in a layer called a photon dominated region (PDR) (see, e.g., van der Werf et al. 1996, 
and references therein). The chemical and radiative processes along with the gas dynamics 
are complex and definitely beyond the scope of this paper. Here, we want to adopt a very 
simple model for such a PDR region and concentrate on the treatment of the continuum 
radiation. Assuming that all kinematic effects of the radiation and of a possible wind are 
influencing only the outer layers of the parent cloud, still the long-wavelength part of the 
radiation is able to reach the dense core deep in the cloud and heat it. The gas atoms, 
molecules, and dust particles convert the UV radiation to radiation at longer wavelength 
which propagates deeper into the cloud and the dense core before being absorbed and re- 
emitted at far- infrared wavelengths (Liseau et al. 1999). For this test application, we will 
neglect the possible influence of a shock front that is formed by the radiation ram pressure 
and serves as an additional source of energy. We also neglect heating from molecular fines 
other than PAH emission. The density and temperature of the dense core is discretized 
on a 3D grid with N — 64000 grid cells. To calculate the incoming energy density of the 
stellar radiation, A'^ traces have to be performed to derive the column density from the star 
to the cell. Then, in each cell, the stellar intensity of the cells can be determined for all 
wavelengths. 

The PDR is discretized by a layer with Np^^ x Np^ji x 1 cells at a distance of 0.5 
pc covering a cloud surface of 0.25 pc^. The emission of dust particles can be treated 
as black body radiation for particles with a size of a > 0.1 /im, making it possible to 
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calculate the energy balance for the dust particles in each cell quickly. Smaller dust particles 
and PAHs absorb and re-emit the radiation time-dependently, substantially increasing a 
numerical treatment. For the purpose of this test calculation, we will therefore assume that 
the conversion factor of stellar energy from the UV to optical wavelengths by the PDR is 
a known value rj and parameterize the PDR spectrum by a power-law in wavelength of the 
form XFx ~ A~° *^ between 5 /im and 18 /xm (Zucconi et al. 2001). To calculate the density 
of the radiative energy absorbed in the cells, iVj,^^ x N traces have to be performed (we 
have chosen NpoR — 40). The interstellar radiation field will mainly contribute to the 
heating in the far-infrared (FIR). Assuming it to illuminate the configuration isotropically, 
we have to perform traces covering all directions. As direction grid on the unit sphere, we 
use an equally-spaced distribution of Nd = 200 points optimized for integration purposes 
using simulated annealing (Steinacker et al. 1996). The spectrum proposed by Zucconi et 
al. (2001) was used for the interstellar radiation field. In total, N^, x N traces have to be 
performed to determine the heating by the interstellar radiation field. The density at any 
location was calculated from the grid using Ist-ordcr interpolation. 

Finally, to obtain a self-consistent temperature distribution, the cell-cell illumination 
has to be considered incorporating A^^/2 — N x Ni traces for the column density between 
the cells, where is the number of iteration to reach convergence of the cell temperatures. 
For cores, the self-heating is restricted to the FIR- and mm-wavelengths range as the dust 
particles have low temperatures around 20 K. This contribution to the cell temperature is 
low, though, since the absorption and re-emission coefficient k"^*(A,x) decreases as in 
this region. Thus, the temperature converges to a solution after only a few iterations. The 
total amount of traces for this example adds up to 2.2 x 10^ encountering integrated optical 
depth up to 10"^. We have used the optical dust properties of Silicate particles of size 0.1 //m 
proposed in Draine & Lee (1984) in this calculation. 

Fig. 3 shows the optical depth t(10 /im) for rays crossing the core parallel to the x-axis 
(y and z are given in arbitrary coordinates in the plane of sky). The optical depth reaches 
values around 500 in the densest parts. The resulting temperature structure is visualized in 
Fig. 4 for the case of external heating by the interstellar radiation field. The contour plot 
shows a cut through the temperature cube in a plane parallel to the plane of sky including 
the point with the temperature minimum. In massive star formation regions, this is the case 
only when no massive star has formed in the vicinity yet. Fig. 5 gives the temperature in 
the same plane when a nearby B6 star and the corresponding PDR is present. While the 
isotropic interstellar illumination yields a temperature structure that is similar to the density 
structure, the beamed illumination by the star and the PDR causes shadowing effects with 
low temperature behind density maxima. This effect will be important for the collapse of the 
core due to self-gravity as the local instability criterion depends on the temperature. The 
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influence of external heating by nearby stars will be investigated further in a forthcoming 
paper. 

5. Summary 

There is a variety of astrophysical applications requiring the calculation of the radiation 
field in complex structures with regions of high optical depth. Current radiative transfer 
codes suffer from the high computational cost to determine the radiation field for such cases 
as the existing solution algorithms are often not suited for high optical depths. In this 
publication, we have suggested a ray-tracing scheme designed for the analysis of complex 
and dense structures based on 3D radiative transfer. Analyzing accuracy, numerical effort, 
and stepsize for the example of continuum radiation, we have found it to be well-suited for 
applications including regions of arbitrary high optical depth and a vast number of required 
ray-traces. For the test case being investigated, the tracer was about 870 times faster than the 
full 5th-order Runge-Kutta tracer with a loss of only one order of magnitude in accuracy. 
We have applied the scheme to a massive molecular cloud core illuminated by a nearby 
massive star. Using a simple model for a photon-dominated region at the surface of the 
parent molecular cloud, we have compared the 3D temperature structure within the core for 
two cases: (i) an illumination by the interstellar radiation field and (ii) an illumination by 
a nearby massive star and a corresponding photon-dominated region. Stellar illumination 
produces shadowed regions of colder gas and dust which will influence the evolution of the 
core possibly forming stars. 

J.S. thanks the Institute for Pure and Applied Mathematics, UCLA, and the Observa- 
toire de Bordeaux, L3AB, for support and beneficial discussions during visits in 2005. 
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Table 1: Typical peak optical depths (A = 550 nm) for dense dusty environments 
environment r 
low-mass molecular cloud core^ 100 
massive molecular cloud core*^ 10^ 
circumstellar disk around a young star'^ 10^ 
dust torus around an active galactic nucleus*^ 100 



"derived from the optical depth at 7 /xm found in Bacmann et al. (2000) 
''Wu et al. (2005) 

"^not measurable yet, standard accretion disk model and common interstellar dust grains are assumed 
''see, e.g., Watabe & Umemura (2005) 
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Fig. 1. — Overview of the physical set-up and solutions for the test case calculation. Top 
left: Optical depth along the ray. Bottom left: Temperature along the ray. Top right: 
Intensity along the ray using a 5th-order Runge-Kutta solver for the ray equation. Bottom 
right: Normalized difference between the solutions obtained from the accurate solver and the 
relative intensity ray-tracer (thin line) and interpolation error (thick line). Regions where 
the solver uses the optically thick approximation are enclosed by vertical dotted lines (z=0.9 
to z=3.2, and around z=6). 
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Fig. 2. — Top left: Relation between the run time and the maximal error between relative 
solution and accurate solution. Top right: Run time as a function of the total optical depth 
for a ray crossing the region. Bottom left: Run time as a function of the density gradient 
maximum. Bottom right: Run time as a function of the maximal gradient in the source 
function. All quantities aside from the maximal error are normalized to the test case shown 
in Fig. 1. 



Fig. 3. — Optical depth of a massive molecular cloud core with a density structure adopted 
from a multi-wavelength 3D radiative transfer modeling of the low-mass core p Oph D at 
A = 10 pm. along the line of sight (x-axis). The local plane-of-sky-coordinates are arbitrary. 
The optical depth reaches values around 500. 
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Fig. 4. — Iso-temperature distribution within the cloud core in a plane parallel to the plane of 
sky that includes the temperature minimum. The core is heated by the interstellar radiation 
field only. The numbers give the temperature in K with a 1 K stepping in the contours. The 
local plane-of-sky-coordinates are arbitrary. 
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Fig. 5. — Iso-tcmperature distribution within the cloud core in a plane defined in Fig. 4. 
The core is heated by a nearby B6 star, a corresponding PDR, and the interstellar radiation 
field. The numbers give the temperature in K with a 1 K stepping in the contours. The local 
plane-of-sky-coordinates are arbitrary. In this orientation, the star and PDR are located 
above the picture. 



